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ABSTRACT 

Following a suggestion that a directed relativistic explosion may have a universal intermediate 
asymptotic, we derive a self-similar solution for an ultra-relativistic jetted blast wave. The solution 
involves three distinct regions: an approximately paraboloid head where the Lorentz factor 7 exceeds 
^ 1/2 of its maximal, nose value; a geometrically self-similar, expanding envelope slightly narrower 
than a paraboloid; and an axial core in which the (cylindrically, henceforth) radial flow u converges 
inward toward the axis. Most 80%) of the energy lies well beyond the leading, head region. Here, 
a radial cross section shows a maximal 7 (separating the core and the envelope), a sign reversal in 
R, and a minimal 7 , at respectively ^ 1/6, ~ 1/4, and ^ 3/4 of the shock radius. The solution is 
apparently unique, and approximately agrees with previous simulations, of different initial conditions, 
that resolved the head. This suggests that unlike a spherical relativistic blast wave, our solution is an 
attractor, and may thus describe directed blast waves such as in the external shock phase of a 7 -ray 
burst. 

Subject headings: hydrodynamics — jets — relativistic processes — gamma-ray burst: general 


1. INTRODUCTION 

Relativistic jets give rise to some of the most spectac¬ 
ular astronomical sources, including 7 -ray bursts (GRB) 
with Lorentz factors T ^ 10^-10^, active galactic nu¬ 
clei (AGN) with T ^ 10, and microquasars with T ^ 
a few. Such jets are also likely to form in systems cur¬ 
rently invisible to us, such as failed supernova explosions. 
Yet, there is considerable uncertainty regarding the ori¬ 
gin, structure, and evolution of relativistic jets. 

One might expect that in the limit of an evolved, ultra- 
relativistic jetted blast wave (henceforth jet), propagat¬ 
ing into a sufficiently weakly magnetized, homogeneous 
medium, the lack of a characteristic length scale would 
lead the structure of the jet to approach some self-similar 
attractor. Such a solution may provide a framework for 
studying jets, and be useful as an approximate descrip¬ 
tion of directed blast waves in astronomical systems, such 
as in the external shoc ks stage of a GR B. 

It has been argued (jGrnzinovI I2QQ7L henceforth GOT) 
that a directed explosion, in which the total momentum 
n is comparable (when multiplied by the speed of light 
c) to the total energy Stot^ such that Stot — He ^ Stot^ 
may have a well-defined self-similar attractor. This dif¬ 
fers from a non-directed {Stot ~ nc < Stot) explosion, 
which is not in full causal contact when highly rela¬ 
tivistic, and thus has no universal intermediate asymp- 
totic. Indeed, non-d ir ected explosions asymptote to the 
iBlandford fc McKed (|1976[) self-similar solution in the 
ultra-relativistic phase only if t hey are initially spheri¬ 
cally symmetric (|Grnzinovll2QQQ[) . 

The struc ture of a relati vistic jet was described qual¬ 
itatively in iRhoadsI (|1999[ ). and quantitatively in GOT, 
where numerical simulations of various initial conditions 
were argued to approach a unique attractor solution. The 
reported jet has a nearly paraboloidal shock front, and, 
with increasing distance from the shock, shows a mono¬ 
tonic decline in the local Lorentz factor 7 , the proper 

ukeshet@bgu.ac.il 


pressure p, and the (cylindrically, henceforth) radial ve¬ 
locity u. 

While these results are promising, the universality, 
uniqueness, and full structure of the attractor were until 
now unclear. The jet structure was derived from simula¬ 
tions rather than by directly solving the flow equations. 
This structure was reported only for the head region, 
within a distance 1 < ^ < 5 of the light signal preceding 
the jet, normalized such that ^ = 1 is the nose position 
(see Eq. ([T8|) for a precise definition of ^). The simula¬ 
tions were reported to approach the attractor very slowly, 
suggesting that they have not fully converged upon it. A 
possible sign reversal in u was reported far downstream, 
probably in the non-converged region, suggesting some 
deviation from a parabolic profile (GOT). 

Several questions have thus remained open. What 
are the full equations describing the self-similar, ultra- 
relativistic jet? These equations were only partly derived 
in GOT. What are the solution or solutions to these equa¬ 
tions? What is the corresponding normalization of the 
self-similarity scaling? Is the solution unique, and if so - 
what provides closure to the equation system? What is 
the structure of the jet, in particular far from the nose? 
Does it show the flow monotonicity reported at the head? 
Is most of the energy contained in the head, or does it 
lie beyond it, in regions not yet reported, and probably 
not converged? 

Here we address these questions. First, we derive the 
equations for an ultra-relativistic, self-similar jet in ^ 
following and supplementing GOT. Analytic and semi- 
analytic constraints on the solution are derived in 0 
In 0 we finally solve the equations numerically, and 
show that their regular solution is unique and satisfies the 
constraints obtained in 0 The results are summarized 
and discussed in ^ The observational implications and 
stability of the solution are deferred to a forthcoming 
publication. 

We assume ///an axisymmetric relativistic flow ex¬ 
panding into a homogeneous medium in fiat space; 
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(ii) that the effects of cooling, electromagnetic fields, and 
self-gravity are negligible; and (in) self-similarity in time. 
We normalize the speed of light, c = 1. 

2. SELF-SIMILAR ULTRA-RELATIVISTIC EQUATIONS 

2.1. Hydrodynamic equations 

Following GO7 and using similar notations, we con¬ 
sider a directed fiow propagating along the positive 2 ; 
direction, where = (t, x, z) are cartesian coordi¬ 
nates in the upstream frame. The fiow downstream of 
the shock is parameterized by the squared Lorentz fac¬ 
tor g = 7 ^, the proper pressure p, and the cylindrically 
radial (Le. perpendicular to the axis of symmetry) ve¬ 
locity u. We thus analyze these downstream properties 
using upstream frame coordinates, henceforth omitting 
the terms upstream and downstream in this context. 

The fiow is governed by the relativistic hydrodynamic 
equation 

= 0 , ( 1 ) 

where = (Au^u^ — g^^)p is the energy-momentum 
tensor, is the four-velocity, is the Minkowski 
metric, and we assumed a relativistic fiuid with an en¬ 
thalpy density of 4p. Assuming axial symmetry in cylin¬ 
drical coordinates (t,r, 0 , z), the four-velocity becomes 
= 7 ( 1 , + The t, z and r components of Eq. ([T]) 

are then 

dt [(4g - l)p] + Adz{qpv) + Ar~^dr{rqpu) = 0 , (2) 

4dt{qpv) -h dz [{Aqv^ -h l)p] + Ar~^dr{rqpuv) = 0 , (3) 
and 

Adt{qpu) -L Adz{qpuv) + Ar~^dr{rqpu‘^) + drP = 0 . (4) 

In the upstream medium (denoted by a tilde), pres¬ 
sure is negligible, so the only non-zero component in its 
energy-momentum tensor is e = Parameter¬ 

izing the shock (subscript 5 , henceforth) surface as the 


zero isosurface of a scalar function 

S{x^) = z - Zs{t,r) = 0 , (5) 

the requirement of continuous energy-momentum fiuxes 
across the shock can be written as 

. (6) 

The t, ^ and r components of these constraints become 
AqsPs{dtZs - Vs A- UsdrZs) - pdfZs = edfZs , (7) 

'^QsVs(^dt^s Vs Usd-pZs^ 1 — 0 , ( 8 ) 

and 

AqsUs(^dfZs Vs H“ Usd^Zs^ fi dj^Zs — 0 . (9) 


The fiow is given by a solution to equations (EHll), with 
shock boundary conditions (IMI)- Notice that in Eq. m 
we assumed that 2 :^ is a function of r^, and not vice 
versa, such that the jet is infinitely long, and monotoni- 
cally widens with z. This excludes possible scenarios in 
which the jet widens near the head but narrows down 
farther downstream, as found in some si mulations {e.g., 
iGranot et al.l I 2 QQQI: iGannizzo et al.l l2QQ4f ) . The assump¬ 
tion can be avoided by parameterizing the shock using 
rs(t, 2 ;), instead of see ^ 


2 . 2 . Ultra-relativistic limit 

In the ultra-relativistic limit, q~^ can be used as a 
small expansion parameter, giving 

v= y/l-q-^ -u^ = 1 - ^ + O (g~^) , (10) 

where we assumed the scaling qu‘^ = 0 ( 1 ), as confirmed 
by the results (see G07 and Q. We define 

wi = 1 qu^ and 1^2 = 1 + 2qu^ , ( 11 ) 

for future convenience. 

It is useful to replace 2 ; by a coordinate C = t — z^ mea¬ 
suring the distance to a plane perpendicular to the jet 
and moving ahead of it at the speed of light. Equations 
(|2]-[4]) then become, to leading order in 

49t(gp) + d(^{w 2 p) + {A/r)dr{rqpu) = 0 , ( 12 ) 

dt{w 2 p) + d(^{wlp/q) + {2/r)dr{rwipu) = 0 , (13) 

and 

AqpdfU + 2pwid(^u + ud(^p + drP + AqpudrU = 0 . (14) 

Note that we added the last term in Eq. m to correct 
equation (18) of G07. 

In a similar fashion, the shock boundary conditions (13- 
| 2 ]) at C = (^s{t,r) = t — Zg become, to leading order in 

" AdtCs + 2{drCs? ’ • 

(15) 

2.3. Self-similarity 

The ultra-relativistic Eqs. dUHia) remain invariant if 
the vector 

V= (Ct,r^C^g,CVCVw^) (16) 

is multiplied by an arbitrary constant, and so they admit 
a self-similar scaling. In D spatial dimensions, the total 
energy of a self-similar solution would scale as 

Stot OC OC C* ^ , (17) 

where a star denotes a typical value at time t. This en¬ 
ergy is constant during the self-similar phase, imposing 
scalings such as C* (x and r* (x . 

These scalings hold for D 7 ^ 3, but in three spatial dimen¬ 
sions they degenerate into exponential functions. Indeed, 
the equation system remains invariant under rescaling of 
the vector V in Eq. (USD, even if the first, t-dependent 
component is eliminated from V. 

The resulting exponential scaling motivates the dimen¬ 
sionless, self-similar, respectively axial and radial coor¬ 
dinates 

f = A^— and 7 = A— , (18) 

r T 

and the self-similar functions describing respectively the 
Lorentz factor squared, the pressure, and the radial ve¬ 
locity, 

Q = h~^q , P = A“^^ , and U = Au . (19) 
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Here, the temporal scaling factor 

t—tf 

A = e—^ (20) 

is large in the self-similar regime, 

( 21 ) 

is the e-fold expansion time, C is a dimensionless con¬ 
stant, and we defined a final time tf when the flow is no 
longer relativistic. The approximation holds for g ^ 1 
and ^ 1 , so (t — t/) is assumed negative and large with 
respect to r. 

The total energy in the jet is related by Stot = T^eEtot 
to the total self-similar energy 

/ nOO nOO 

EdV = Sir T]dri d^ QP , ( 22 ) 

Jo J^siv) 

where the self-similar energy density is ^ = 4:QP. This 
fixes the constant C = . 

Rewriting Eqs. (inHia) using these definitions finally 
gives the t, 2 ) and r components of Eq. m in the self¬ 
similar regime. 



+ (23) 

(2 + a*)(«ijP) - dfiwlP/Q) - = 0 , (24) 

and 

(Udr, + 1 - d^)U +^d^U+ P = 0, (25) 


where we defined V’ = log(^^/^r?) for brevity. The bound¬ 
ary conditions are written in self-similar form on the sur¬ 
face ^ = ^s( 7 ?) = (A^/t)Cs, as 


Qs = 


86 


1 

+ 26" ’ 



and f/g = 6 • 


(26) 

The scaling m implies that ici = 1 + QU^ and W 2 = 
1 + 2QU^ can be written in a self-similar fashion, in a 
form identical to Eq. m- Note that we added the first 
term in the parenthesis of Eq. (ESI) to correct equation 
27 in GOT. 


2.4. Equation properties 

The self-similar Eqs. (|23H25|) and boundary conditions 
(ESI) remain invariant if A is multiplied by an arbitrary 
constant. This constant can in principle be chosen as 
(— 1 ), indicating that Q and P are symmetric, while U is 
antisymmetric, under the reflection r] (—(g) across the 
axis of symmetry (henceforth, the axis). 

Without loss of generality, henceforth we choose this 
arbitrary constant such that the nose, i.e. the very head 
(subscript h, henceforth) of the jet is at 

^ = 6 = Uv = 0 ) = 1 . (27) 

Equation ([26|) then implies that at the head, the flow 
parameters are given by 

Qh = ^ , = ^ ^ Uh=0 , and Eh = ^ . (28) 


Here and below, we define Fg, Fa and Eh as the field 
E G {Q, P, t/, E} evaluated at the shock, on the axis, 
and at the head, respectively. 

Notice that Eqs. (j23]-[25j) are, in addition, unchanged 
if P is multiplied by an arbitrary constant, but this con¬ 
stant is fixed by the boundary conditions (|26]) . 

The equation system involves three first-order partial 
differential equations, along with three boundary con¬ 
ditions, for the three fields {Q^P^U} living in the two 
dimensional ^-77 space. In addition, the equation system 
depends on the unknown one-dimensional function ^s(g), 
so one additional constraint may appear to be missing. 
As we show below, the system is closed, and the shock 
profile ^s{v) is fixed, by requiring that the solution be 
regular, even without specifying boundary conditions far 
downstream. 

The oriented derivative along the shock may be defined 
as 

ds = drj + C{v)di ■ (29) 

Eor a shock profile that is monotonic, in the sense 
that the function ^s{v) monotonically increases, one 
may alternatively write the derivative in Eq. (|2^ and 
the boundary conditions (| 2 fi|l using the parametrization 
? 7 s(^), instead of Cs(^)- This more general parametriza¬ 
tion is used in 0 

If the shock profile is known, or is postulated, one may 
use Eq. ([29|) to turn the self-similar Eqs. (|23H26|) into 
an infinite series of ordinary differential equations for in¬ 
creasingly high-order field derivatives, in either the ^ or 
the 77 direction. This is utilized below, in 0and in m 

Eor a monotonic shock profile, C'(g) > 0, Eq. (j26j) 
implies that Us > 0 , i.e. near the shock, the radial flow 
is always directed outwards. Combining Eqs. (|25]) . (|26]) 
and (EH) indicates that 

d^Us = 0 , (30) 

which may be interpreted as a vanishing self-similar ra¬ 
dial acceleration. Note that the term we added to equa¬ 
tion 27 of GOT, in order to obtain Eq. (|25]) . is essential for 
recovering this property. In particular it implies, using 
Eq. ([29]), that dr^Us = Cs • 


3. SEMI-ANALYTIC CONSTRAINTS 
3.1. Overview of the solution 

Expand the shock profile about the nose of the jet, 
fs = 1 + ^ 2 ^^ + 0{rf)^ where ^2 and ^4 are numer¬ 

ical constants, and we used reflection symmetry to omit 
odd powers of 77 . The coefficient ^2 is unlikely to vanish, 
as confirmed numerically below, so the head of the jet 
is approximately a paraboloid. Define the logarithmic 
shock profile slope. 


^ din fs 
din 77 


(31) 


such that at the nose (3^2. 

A parabolic, (3 = 2 shock profile is stationary, in the 
sense that a point {r, Cs(^) oc is mapped at a time At 
later onto {Ar, A^Cs (x (Ar)^}, with A = e^^E^ in non 
self-similar coordinates, a self-similar (3 < 2 {i.e. wide 
jet) profile narrows down in time, eventually approach¬ 
ing (^s ^ whereas a /3 > 2 (narrow) shock gradually 
widens toward (s ^ r‘^. 
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Far from the head, the self-similar shock profile must 
become narrower than a paraboloid, i.e. /d > 2, because 
/3 < 2 leads to nonphysical results, in particular a diver¬ 
gence of the energy in Eq. (j22j). Moreover, P < 2 leads 
to only non-real valued solutions (see while P = 2 
leads to a divergence of Q (see ^3.31 and Figured]), and 
cannot be matched to the axial region (see ^Bj). 

At large distances from the head, one may approximate 
/3 as a constant, up to logarithmic corrections. For such 
power-law behavior, oc 77^, an additional, geometrical 
self-similarity (GSS) in the plane may emerge far 
from the head and from the axis, with the self similar 
parameter (GOT) 

Xo = , (32) 

as discussed in ^3.21 

A GSS scaling in the relevant, P > 2 regime (demon¬ 
strated in ^3.41 and in Figure [2|) breaks down near the 
head and near the axis, as equation terms that do not 
follow the GSS scaling become large and can no longer 
be neglected. Denote the surface along which the axial 
terms break the GSS scaling as ^ = Cc(^)- One may de¬ 
fine this surface as the boundary between an axial, or 
core, region at ^ and the GSS region at ^ < <^c- In 
0 we show numerically that this boundary can be asso¬ 
ciated with a maximal value of Q{r])^ giving — 20r]^ . 

In our P > 2 regime, Eqs. (jSSHsOj) imply that drills = 
> 0, so near the shock the (positive) radial velocity 
increases radially outwards. At large radii, P > 2 implies 
that Q and P similarly become monotonically larger as 
one approaches the shock, such that d^{Qs^Ps} < 0 and 
dr]{Qsj Ps} > O5 ns seen by solving Eqs. (|234B6|) . with 
the aid of Eq. (1221, for these derivatives at the shock. 

This behavior, in which F G {Q^ E} > 0 increases 
monotonically toward the shock, namely 

d^F < 0 and > 0 , (33) 

is henceforth referred to as flow monotonicity. The sim¬ 
ulations of GOT show this type of behavior throughout 
the reported, 1 < ^ < 5 regime. It is therefore natural 
to ask if flow monotonicity persists throughout the jet. 

Within the GSS regime, P > 2 solutions are found to 
transition at some Xo = Xu, from U (xo < Xu) >0 t owar d 
the shock, to I/(xo > Xu) <0 toward the axis (see ^3.4|) . 
Indeed, far from the head, Qa{^'^ ^h) ^ approaches 
a power-law along the axis, with a < 1, which, as we 
show, implies that dr^U{r] = 0) ^ (a — 1) asymptotes 
to a ne gativ e constant (except in the special case a = 
1; see gSIH). This is consistent with a 1/ < 0 inflow 
emanating at Xu within the GSS regime, and extending 
all the way to the axis. Thus, the radial velocity reverses 
sign, corresponding to a radial inflow near the axis, and 
a radial outflow near the shock. As the radial velocity 
vanishes along the axis, Ua = 0, this indicates that tf 
cannot be monotonic. 


To see that Q cannot be monotonic either, recall that 
along the axis, Qa cx: declines with increasing ^ no 

faster than In contrast, along the shock, Eq. (j26|) 
with p > 2 implies that Qs oc ^-2(1-^ ) declines faster 

than . Hence, there must be a region where Q declines 
as T] increases toward the shock, ruling out monotonicity 
in Q. This behavior is confirmed in the GSS regime (see 
^ 3.31 and ^ 3 . 4 [) . 

Nevertheless, flow monotonicity does manifest in the 
head region. Requiring such monotonicity near the head 
implies that the jet cannot be too narrow, constrains the 
flow in the head region, and limits the extent of this 
region to ^ 5 (see ^3.61 and Figure [3|). 

Although the shock profile is not directly constrained 
by the system of equations, and no far downstream 
boundary condition is imposed, we find that only a 
unique profile avoids a divergence of the flow. For ex¬ 
ample, if the jet is too wide, g = 7^ becomes negative on 
the axis (see m, which is both nonphysical and leads to 
divergencies. Indeed, as shown in 0 numerically solv¬ 
ing the equations and requiring a regular jet picks out a 
unique flow solution (see Figures HH^), which agrees with 
all the features mentioned above and derived quantita¬ 
tively below. 

The above arguments indicate that a regular self¬ 
similar jet has a unique solution, composed of three dis¬ 
tinct flow regions: a monotonic head region (^ ^ C^), an 
effectively one-dimensional, GSS envelope (C^ ^ ^ Cc), 

and an axial, or core, region (^ ^ Cc)- A radial inflow 
encompasses the core and the inner envelope regimes. 
Only the head and the outer part of the envelope, which 
harbor a radial outflow, show a monotonic flow. 


3.2. Geometric self-similarity 


At large distances from the head of the jet, the shock 
profile may be approximated by a power law, in which 
case the temporally self-similar Eqs. (|23H26|) may be fur¬ 
ther simplified. Using a scaling parameter xo such as 
defined in Eq. I|32|). these equations can be cast in a ge¬ 
ometrically, and not only temporally, self-similar form. 

It is useful to introduce a slightly different GSS param¬ 
eter. 


X = 


At]^ 


(34) 


where the unknown normalization constant A is intro¬ 
duced such that the position of the shock in the GSS 
regime can be taken as x = Xs 5 with a constant Xs which 
we choose as Xs = 1- 

The case P <2 does not lead to physical GSS solutions, 
as mentioned above and discussed in ^ Therefore, here 
we focus on P >2. 

The flow equations can then be written approximately 
as functions of xmI the parameters are rescaled as 


Q(x) = , P{x) = ^^P , and U{x) = A . 

Equations (|23H25|) now become, respectively, 

“ 2(/3-1)^2 + Az^iixQ'Py + SQV] = 0 


(35) 


{W2'P)' 


4 


- 4)QVU - PxiQPUy 


(36) 
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(^) “ 2/3xKPW)' - 2(/3 - 2)w^VU - ^2 = 0 , (37) 


and 


(^1 - 2^xQU)U' - (2 - ^2)(/3 - 1) + (^/ - /3x) 


4^ 


7^' 


2p ^2(/3-1)^2 




2(/3-2)x^Q 


(38) 


where the scaling (j35j) maintains rci = 1 + QZ//^ and 
1^2 = 1 + 2QU‘^ in GSS form. Note that in Eq. (|38|) . we 
added the second parenthesis, and the second term in the 
first parenthesis, to correct equation 36 of GOT. 

The last two terms in each of the above three equations 
(|36H38|) do not follow the GSS scaling. However, the first 
of the two terms is negligible far from the head, and the 
second is negligible far from the axis. Note that the last 
term in each equation vanishes in the special case [3 = 2. 
GSS behavior is therefore expected to emerge far from 
the axis (for p > 2), or even on the axis but far from 
the head (for /3 = 2). Far from the axis (and thus also 
from the head), the shock boundary conditions on y = 1 
asymptote to the GSS form 

2s = ^ , = ^ , and Us=l3 . (39) 

Finite energy solutions exist only for /3 > 2. Indeed, 
in 21 we numerically find that the full (self-similar but 
generally non-GSS) solution asymptotes to j3 2.02 far 
from the axis. Before addressing such j3 > 2 solutions, we 
first discuss the special case /3 = 2, for which an analytic 
flow solution can be found. Although the GSS equations 
are in principal valid in this case even as ^ 0, the 
solution itself is shown to diverge near the axis. 


3.3. Analytic GSS solution for jS = 2 

As mentioned above, in the special case /3 = 2, the last 
term in each of Eqs. (|M|) - (|38|| vanishes. The equations 
can then be solved analytically far from the head, where 
the terms are negligible. For the boundary conditions 
ii, the solution is 


u =6x - 51 - A»rX(4x2 + 27) + 12x2 - 54 - , 


Q =(8x7/ - 27/2) 



4 

7/-4X 


W \ 

n-u) 


dx 


(40) 


where for brevity we defined gi = [4%^+ 3(^2— 3)^ 
and g 2 = [27 + 32%^ -h 8%^ (27 + 16%^)^/^]^/^. This solu¬ 
tion is shown in Figure [11 as an azimuthal cross section 
through the jet. 

Far from the shock front (y ^ 1), the leading terms in 
this solution are 


“ 27 64x2 


V- 


and 


0.046 
X 

8X 


^ 21rfA'^ ’ 

0.046 

^ 277 / 3^2 


(41) 


Hence, although Eq. (gni) provides an exact solution to 


the GSS equations when ^ 0, and so gives an asymp¬ 
totic solution for f, ^ this solution is not physical on 
the axis, where Q diverges. This local divergence is more 
severe than, and is not directly related to, the global log¬ 
arithmic divergence of the total energy, associated with 
the marginally wide (/3 = 2) shock profile. 

Nevertheless, the (3 = 2 solution does provide an ade¬ 
quate analytic approximation of the jet in the outer en¬ 
velope region, where y is not too large. For example, 
in a radial cross section, it shows a minimal Q where 
0 = dr^Q (X PxQ' + 2(/3 — 1)Q. This occurs at y 1.61, 
or equivalently at a fraction 

f = ^ = x“'/^ 0.79 (42) 

Vs 

of the shock radius, close to the numerical value found 
for the full solution in 21 Moreover, its P and U profiles 
(see Figure [T]) agree qualitatively with the full solution 
even in the head (Figures [3] and |S|) region, as well as with 
the simulated (GOT) head. 

3.4. GSS with 13 = 2.02 

As the 13 <2 GSS regime is ruled out in ^3.3l and in ^Al 
here we consider p slightly larger than 2. The numerical 
solution in 2liii4icates that far from the head, the shock 
indeed approaches a /3 2.02 profile. Accordingly, we 

now derive the GSS solution for [3 = 2.02, but point 
out that the qualitative features of the solution are not 
sensitive to the precise value of [3 in the 2 < /3 < 3 range. 

In this range of (3, the GSS Eqs. (|364[39|) form a closed 
system of ordinary differential equations, which can only 
be solved numerically. The solution for /3 = 2.02 is shown 
in Figure |2j 

As the figure shows, while the P profile is not qualita¬ 
tively changed with respect to the analytic /3 = 2 case, 
the profiles of Q and [/, and subsequently also of F, are 
significantly altered. 

Most importantly, Q no longer diverges near the axis, 
and never becomes a function of g alone. Such /3 > 2 
solutions are therefore physical, unlike the diverging /3 = 
2 profile, and can in principle be matched to the near-axis 
solution. 

Interestingly, the U profile is not everywhere positive, 
i.e. the radial flow is not everywhere an outflow. Unlike 
the (3 = 2 case, U becomes negative at large y, corre¬ 
sponding to a radial inflow converging on the axis. For 
/3 = 2.02, the transition occurs at y 20.3, or equiva¬ 
lently at / 0.23. The exact location of the transition is 

sensitive to the precise value of (3. For example, P = 2.04 
gives / 0.32. 

The minimum of Q along a radial cross section is found 
at / 0.77, not far from the corresponding minimum 

Eq. 221) in the p = 2 case. This result is less sensitive 
than the U = 0 contour to the precise value of /3, giving 
for example a similar, / = 0.75 value for p = 2.04. 
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Fig. 1.— Self-similar jet in the GSS regime, for the special, analytically solvable case /3 = 2. The azimuthal cross section shows the self- 
similarly sc aled Lorentz factor squared Q, proper pressure P, radial velocity P, and energy density E (see labels). Colormaps (cubehelix; 
IGreen|[2011|j and contour s (in interval factors of 2 for Q, P and P, and in factors of 4 for E) show each quantity, scaled by the shock width 
normalization A (see Eq. 1341) . with P, Q and E normalized also by their head values. Note the divergence of Q (and so, also of E) toward 
the axis. 


3.5. Axial expansion and its monotonicity constraint 

As shown in ^3.3l and ^3.41 above, the axial region of the 
jet is distinct from the GSS envelope. It is al^ distin¬ 
guishable from the head region, as shown in ^3.61 below. 
Indeed, the axial structure shows that the monotonic na¬ 
ture of the flow near the head, in which the fields F in¬ 
crease monotonically toward the shock in both the (—0 
and T] directions, cannot hold far beyond the head region. 

To see this, expand Q and P near the axis in even 
powers of r^, 

Q{i,v) = QaiO + Q2i0v^ + + ■■■ (43) 

and 

P{C V) = PaiO + P2{iW + ^4(e)y + • • • , (44) 

and expand U in odd powers of 77, 

= UliOv + + ... . (45) 


Here, the Fn are numerical factors. An 0{r]) expansion 
of the flow equations near the axis now shows that 


U,(f) = d,U((,v = 0) (46) 

Ql ^ l + 4CQa 

For large ^ ^ one can approximate the axial be¬ 
havior of Q as a power-law, Qa — where Qao is 

a constant, so Eq. (gSl) gives 




7 _ 6(4a - 3)Qao^ 

2 4(5a0 


. (47) 


If a > 1, then Ui{^ ^ ^h) diverges due to the third 
term. In addition to being nonphysical, this also breaks 
monotonicity, as at the head Ui is finite. If a < 1, then 
Ui{^ 00 ) = (a — 1 < 0. This rules out monotonicity 
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Fig. 2. — GSS profiles for the numerically-motivated case /3 = 2.02. Notations are as in Figure [T] Notice the radial inflow {U < 0) 
emerging near the axis, characteristic of /3 > 2. 


as well, and is consistent with the U < 0 inflow inferred 
near the axis from the GSS analysis (see D. 

A fully monotonic flow near the axis is thus possible 
only for the special case a = 1, such that 


Near the head, the shock profile is nearly parabolic, so 
one may approximate + ^77^. Here, the arbitrary 

constant A coincides with that defined in Eq. (j34|) . for 
P ^ 2. Equations (j23l26l29[) then yield 


1 — 2Qqo 

4QaO + ISQao 


= const. > 0 , 


(48) 


which also requires <1/2. We conclude that only a 
Qa profile with a = 1 and QaO <1/2 can yield a fully 
monotonic behavior near the axis. However, as men¬ 
tioned in 93.11 such monotonicity cannot persist radially 
out to the shock, as this would contradict the faster de¬ 
cline of Qs{C) according to the (3 > 2 shock boundary 
conditions. 


3.6. Monotonic head region 

Near the head, previous simulations (GOT) and our nu¬ 
merical solution (g)) suggest a monotonic behavior. This 
has several interesting implications. 


d^Q c, and d^Q ^ , (49) 

SO imposing monotonicity would imply that the jet can¬ 
not be too narrow, i.e. that A < 2/3. The jet cannot 
be too wide, either; A < 0.1 (A < 0.3) can be shown 
semi-analytically (numerically) to lead Q to vanish near 
the head; see El 

One may attempt to solve the self similar equations 
for a monotonic flow. To do so, we numerically minimize 
the sum of the squares of the left hand sides of Eqs. dn- 
l25)) . while imposing the boundary conditions |26l and, in 
addition, constraining the fields to be monotonic. This 
leads to an approximate solution, illustrated in EigureO 
The resulting monotonic profiles qualitatively resem¬ 
ble the structure of the head found numerically and in 
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Fig. 3.— Head structure, found by approximately solving the self-similar equations under the monotonicity constraint (see JMJ. The 
profiles are in approximate agreement with the GOT simulations (dot-dashed contours). Notations are as in Figure [T] Here, the shock 
profile parameters /3 and A are determined self-consistently by the solution. 


GOT, as shown in the figure. We find such solutions only 
^ ^ 5, beyond which the mo notonic fields diverge. 
This suggests that — 5 roughly marks the edge of the 
monotonic head region. 

An approximate, monotonic description of the flow in 
the head region may be found using the axial expansion 
in Eqs. (11313) ; see m The resulting approximation is 
plotted on top of the numerical solution of the head re¬ 
gion in Figure[ 6 l using an approximate shock profile ^s{v) 
inferred from the numerical solution of 0 the figure 
shows, the expansion fits the results rather well near the 
head. 

4. NUMERICAL SOLUTION 

In order to solve the self-similar flow Eqs. (|23H26|) nu¬ 
merically, we expand the shock profile r]s{^) to increas¬ 
ingly high order, and find the optimal jet solution at each 
order, in two different methods. Following the arguments 
of ^ an optimal solution is defined as the most regular 
solution to the flow equations, i.e. the solution diverging 
farthest either from the head or from the shock. The 
results of this procedure indicate that a unique solution, 
which remains regular infinitely far from the head and 
from the shock, indeed exists. 

4.1. Method 

First, we map the (azimuthal cross section of the) jet 
onto the unit square, 0 < {p, a} < 1 , through the trans¬ 
formation 

p = 1 — — and cr = 1 — % . (50) 

Vs C 

This maps the shock and the axis respectively onto p = 0 
and p = 1. The head and infinite downstream (^ ^ oo) 
are similarly mapped onto a = 0 and a = 1. 

Next, the shock profile is parameterized as 

6 = + (51) 

where we choose this functional form, using rf instead of 
T] in the parenthesis, in order to obtain better behaved 
functions. We expand the a-priori unknown function 
/3{a) to order n, and switch from a ^s{v) parametriza- 
tion to the more general, Vs{0 description of the shock, 
such that 


Consider the solution for a given order n. The shock 
profile is defined by the (n + 2 ) undetermined parameters 
A and eo, ei,..., Given some choice of these param¬ 
eters, one can integrate the equations in two different 
methods: 

1 . Start from the (a = 0 ) head boundary conditions, 
and advance toward (a = 1 ) downstream infinity. 

2 . Start from the (p = 0 ) shock boundary conditions, 
and advance toward the (p = 1 ) axis. 

For finite n, both methods eventually fail, as the fields 
diverge at some finite cfmax (or equivalently ^max) in 
method 1 , and at some finite pmax (or equivalently Prnin) 
in method 2 . Indeed, the shock profile must be fine-tuned 
in order to delay the divergence, and uncover a larger 
fraction of the jet. 

We use both methods, independently maximizing ^rnax 
and pmax at each order n by scanning the (n + 2 ) dimen¬ 
sional phase space of the shock profile parameters. Thus, 
we identify the best approximation to the shock profile 
at every order. 

The results of an order n = 2 scan are presented in 
Figure m for the maximization of both ^max (left panel) 
and pmax (right panel). In order to project our four¬ 
dimensional scan onto a two-dimensional figure, here we 
set eo = ei = 0, and vary only A and € 2 . It is useful 
to define the parameter P 2 through €2 = (—1 + / 32 / 2 ); 
this corresponds to a paraboloid, ^V^ profile 

for ^ ^ C/i: smoothly transitioning into a GSS-like, — 
(G +rsj profile as ^ ^ oc. 

Both methods show that the parameters A (0.53- 
0.54) and ^2 — (2.02-2.03) provide the most accurate 
shock profile at this order, in the sense that the diver¬ 
gence of the flow takes place farthest from the head or 
from the shock. The full (n + 2) dimensional optimiza¬ 
tion process converges with n on a unique shock profile, 
with similar A and (32 values, as discussed next. 

4.2. Convergence 

Method number 1, in which ^max is maximized, is much 
simpler and faster computationally than the Pmax niaxi- 
mization of method 2. Therefore, after demonstrating 
that both methods give similar results at low orders, 
we pursue high orders using only the ^max maximiza¬ 
tion method 1 . Table [T] summarizes the results of this 
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Fig. 4. — Independent scans of ^max (method 1; left) and pmax (method 2; right) in the A-^2 phase space, for a second-order shock 

profile = (1 -h A? 7 g ED for details). We find the maximal ^max for {A ~ 0.54, /32 — 2.027}, and the maximal pmax 
for {A ~ 0.53, /32 ~ 2.020}. 


method, providing the parameters of the optimal shock 
profile at each order n. 

The odd n terms do not significantly modify the so¬ 
lution or increase the values of ^rnax- We therefore set 
these terms to zero, and use only the even n terms. This 
is equivalent to taking (3{a^) instead of (3{a) in the shock 
parametrization Eq. m- 

For convergence tests, we thus define N = 1 n/2 as 

the effective order of the shock expansion. Our highest 
order, n = 16, corresponds to = 9, and thus involves 
searching for the maximum of ^max in a 10 (including A) 
dimensional parameter space. Due to the high dimen¬ 
sionality, for orders n = 6 and above, our maximal ^max 
values should be considered as lower limits, as noted in 
the table. 

Figure [5] shows a convergence plot of our results with 
the inverse effective order N~^ of the expansion. As 
the table and figure show, the maximal ^rnax (disks in 
the figure) monotonically increases with n, and sug¬ 
gests convergence near our resolution limit (anticipated 
at imax ^ 50). The shock profile is well behaved, in the 
sense that low order terms do not change significantly 
as higher order terms are added. 

As the order n increases, eo ^ 0 , corresponding to the 
expected paraboloid head. The largest coefficient in the 
high order /3(cr) expansion is 62 — 0 . 011 , corresponding 
to p 2 — 2.022. This dominates the deviation from a 
parabolic profile, and was therefore chosen as the term 
scanned in Figure IH 

Far downstream, as a ^ 1, the shock profile ([52]) ap¬ 
proaches a power-law, oc 77 ^°°, suggesting a GSS scal¬ 
ing with 

n 

^ = ^oo = 2 + ^e,- . (53) 

j=0 

The values of Poo found at different orders n are shown 
in the table and in Figure [5] (diamonds). These results 
suggest that the shock profile converges, as n ^ oc, at 
Poo 2.02. 

The energy E{< ^rnax) is shown in Figure [5] (squares 
plotted against left and bottom axes) to converge slowly. 


m 


0.0 0.2 0.4 0.6 0.8 1.0 



0.00^ —^^^^^^^^^^^^^^^^^^L 

0.0 0.2 0.4 0.6 0.8 1.0 


Ar-i=(l+n/2ri 

Fig. 5. — Convergence plot (with left and bottom axes) of the 
shock profiles rjs in Table ^ as a function of the effective in¬ 
verse order N~^. As N (or equivalently n) increases, ^max (blue 
disks; arrows designate lower limits) increases and approaches our 
resolution limit, 3 00 (red diamonds) converges at ~ 2.02, and 
Etot{< ^max) (green squares, shown also with the right axis) ex¬ 
ceeds ~ 1.4. A better estimate of Etot is obtained by plotting 
Etot{< 0 for a high (n = 10) order expansion (solid curve, plotted 
with upper axis), indicating that Etot — 2.1. These asymptotic 
estimates are based on extrapolations (dashed curves) to n —)■ 00 
(lower axis) or ^ ^ 00 (upper axis). 

The figure also shows (solid curve with right and upper 
axes) a better method to estimate the total energy of the 
jet, by extrapolating to ^ ^ 00 the E{< profile of a 
solution with a fixed high n . The extrapolated result of 
this (good) fit is Etot — 2 . 1 ; most (^ 80%) of this energy 
lies well beyond the head region. 

4.3. Solution existence and uniqueness 

Our convergence tests suggest that the procedure out¬ 
lined above can in principle be continued indefinitely, 
to arbitrarily high order, such that a solution in the 
n ^ 00 limit exists. Namely, given sufficient compu¬ 
tational power, the shock profile could be adjusted such 
that the integration be successfully carried out to arbi¬ 
trarily large ^max and pmax^ 
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n 

^max 


eo 

£2 

64 

ee 

es 

eio 

ei2 

ei4 

ei6 


0 

12.14 

0.523T 

0.0323 









2.064T 

2 

15.93 

0.5398 

0.0001 

0.014 








2.0282 

4 

18.84 

0.5489 

-0.004T 

0.0155 

0.002 







2.0255 

6 

>20.11 

0.5433 

0.0006 

0.0113 

0.0036 

-0.0015 






2.02T9 

8 

>22.IT 

0.544T 

-0.0004 

0.0118 

0.0019 

0.0010 

-0.0018 





2.0253 

10 

>24.T9 

0.545T 

-0.0001 

O.OIOT 

0.0016 

0.0012 

-0.0011 

-0.0012 




2.0221 

12 

>24.80 

0.5452 

-0.0001 

0.0110 

0.0015 

0.0011 

-0.0010 

-0.0011 

-0.0001 



2.0226 

14 

>24.98 

0.5454 

-0.0002 

0.0109 

0.0016 

0.0011 

-0.0012 

-0.0010 

0. 

-0.0001 


2.0222 

16 

>25.28 

0.5454 

-0.0001 

0.0109 

0.0016 

0.0012 

-0.0012 

-0.0010 

0. 

-0.0001 

-0.0001 

2.0225 


TABLE 1 

Optimized shock profile (see Eq. [52j for increasing expansion orders n. 


S5- 



R- 



R- 



1.0 


R- 


ElC. 6.— Head region of the numerical jet solution for a high (n = 10) order expansion of the shock profile, as listed in Tabled The 
analytic model (dashed contours; see © for the corresponding power-law shock profile (/3 = 2.02 and A = 0.53) provides a good fit near the 
head. The profiles qualitatively agree with the simulations of GOT (dot-dashed contours) near the head, although the GOT jet is narrower. 
Notations are as in Figure [T] but here contour values are mostly chosen to match those of GOT, for comparison purposes. 
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Fig. 7. — Full numerical solution for a high (n = 10) order expansion of the shock profile. The difference between the shock front (thick 
solid curved) and a paraboloid (thick dashed blue) is small, but noticeable. Far from the head, the minimal Q is found in the envelope at 
/ ~ 0.75 (dotted cyan curve). Curves with non-monotonic Q first appear at ~ 3.8, roughly located at / ~ 0.17 (dot-dashed red curve); 
this may be regarded as the boundary between the core and envelope regions. An inflow, U < 0 region appears at ~ 5.5, confined to 
/ < 0.33 (dot-das hed yellow curve; the narrow U > 0 stripe along the axis is probably a numerical artifact). The dot-dashed curves are 
plotted using Eq. (El, with replaced by or ^u- 

















12 


Keshet & Kogan 


The existence of a solution extending infinitely far 
downstream is further supported by the GSS analysis of 
93.21 - 93.41 in which a solution for the envelope was indeed 
shown to exist out to ^ ^ oo. However, in the absence 
of a matched solution for the core, we cannot rigorously 
prove the existence of the full jet solution much beyond 
our present, ^rnax ^25 limit. 

Even if we assume that the solution exists, there is no 
a-priori guarantee that it is physically meaningful, be¬ 
cause in realistic scenarios the assumption g' ^ 1 of an 
ultra-relativistic flow eventually breaks down far down¬ 
stream. As the downstream flow is subsonic with respect 
to the shock, the solution could in principle be completely 
altered once the unavoidable subsonic transition is in¬ 
corporated. The physical relevance of the solution would 
however be assured, if it is a sufficiently strong attractor. 

Assuming that the solution exists, it does appear to be 
unique, because (i)we do not identify any other solution 
at low order n < 4, where the phase space can be thor¬ 
oughly searched for additional solutions, as illustrated 
in Figure m (^zi^such a search is double-checked using 
the two scanning methods described in 94.11 found to 
agree well with each other; (iii)hi^ dimensional scans, 
although not as complete, show no evidence of other so¬ 
lutions; and (iv)i]ie semi-analytic constraints of ^ limit 
the phase space of possible solutions. 

Numerical simulations {e.g., GOT) provide the best ev¬ 
idence that the solution indeed exists, is unique, is phys¬ 
ically meaningful, and even behaves as a strong attractor 
and is likely to be stable, at least against axisymmetric 
perturbations. For, in such simulations, roughly tuned 
for external GRB shock parameters, the self-similar so¬ 
lution is found to emerge from different initial configura¬ 
tions, involving various relativistically moving blobs. 

Finally, note that although in ^and 0we assumed an 
infinite jet with 2 ;^ = r), the numerical solution here 

was derived using Eq. (©, thus essentially assuming an 
z) shock profile. Hence, we also search for, and 
rule out, self-similar jet (z.e. relativistic, directed blast 
wave) solutions with non-monotonic r]s{C) profiles, such 
as pinched jets. 

4.4. Jet structure 

The numerical solution to the self-similar structure of 
the jet is shown, for the head region < 5) in Figure 
[HI and for the full range available (1 < <^ < ^rnax) in 
Figure El using a high (n = 10) order expansion of the 
shock profile, optimized in method 1. The corresponding 
shock profile parameters are listed in Table [TJ 

As Figure [ 6 ] shows, the numerical results in the head 
region are fairly well fit by the monotonic head model 
derived in 93.61 (dashed curves). It is also in qualitative 
agreement with the GOT simulations (dot-dashed curves), 
although the latter correspond to a somewhat narrower 
shock profile. 

In Figure El the slight deviation of the shock profile 
from a paraboloid shape (thick dashed curve) becomes 
apparent far from the head, as [3 approaches its asymp¬ 
totic value ydoo — 2 . 02 . Indeed, here the envelope of the 
jet qualitatively resembles the non-mo notonic, p = 2.02 
GSS profile of figure O 

More quantitatively, a radial cross section at large ^ 
shows the minimum of Q at / = rj/rjs — 0.T5 (dotted 
curve), and the sign flip oi U at f 0.33 (dot-dashed 


curve in the U panel). These values deviate somewhat 
from those expected for P = 2.02; both would agree with 
(3 = 2.04 (see 9^ . This is not surprising, as our nu¬ 
merical solution only reaches rjmax — 6.5, which may not 
be sufficiently large to show the asymptotic GSS behav¬ 
ior. Moreover, the 1/ = 0 contour emerges from the axis 
only at ^ 5.5; so it may not have converged by 

^max — 25. 

The breakdown of Q monotonicity is already evident 
around ^ 4, and perhaps even closer to the head, 

but the precise location at which the ’shoulder’ in Q {e.g., 
the wiggle in the QIQh = 0.3 contour shown in Figures 
[6] and El emerges from the axis is not well converged. 

At large this feature corresponds to a maximum of Q 
in a radial cross section, located roughly at / = /c ^ O.IT 
(dot-dashed curve in the Q panel of Figure E])- This may 
be regarded as the boundary between the core and enve¬ 
lope regions. Note that, although the transition appears 
to take place at a constant /, this did not have to be the 
case, as the GSS scaling does not apply in the core. 

The above results support the partition of the jet into 
three distinct regions, as anticipated in 0 a head re¬ 
gion for ^ a core region for ^ ^ Cc and a GSS 

envelope for ^ > Cc- The boundary of the head region 
corresponds to the breakdown of monotonicity, 5, 

and is related to the onset of non-monotonic Q and neg¬ 
ative U behavior, at and ^t/. The boundary between 
the core and the envelope can be identified as fc — 1/6, 
corresponding to Xc — 36 or equivalently — 2677 ^. 

5. SUMMARY AND DISCUSSION 

The equations governing the structure of a self¬ 
similar, ultra-relativistic, directed blast wave are derived 
(Eqs. [23H26|) . following and correcting GOT. A numeri¬ 
cal analysis (21) suggests a converging (Figure [5]), unique 
(Figure 11]) jet solution (Figures [6] and E]), which qualita¬ 
tively agrees with previous simulations (GOT) of the head 
region, and with our semi-analytic study (^. 

The jet can be broadly partitioned into three distinct 
regimes: a head region (^ < ^^ :^ 5), an axial core {^g < 
^ — 2077 ^), and an envelope (^ > ^c)- The highest 

Lorentz factors, 7 > 7 /i/ 2 , are found in the head, most 
of the energy lies in the envelope, and the core contains 
an axial inflow originating from the envelope. 

In the head region, Q > (Q/i/4), P > {Ph/lO), U > 
0 , the shock profile is very close to a paraboloid, and 
the flow is monotonic, in the sense that Q, P, and U 
monotonically decrease away from the shock in both the 
^ and the (— 77 ) directions. This region, analyzed in 93.61 
and in is strongly constrained by monotonicity {e.g., 
Figure [ 3 ]); it qualitatively agrees with the head structure 
reported based on simul ation s in GOT. 

The core, analyzed in 93.51 shows a monotonic behav¬ 
ior in Q and P. The radial velocity U is negative here, 
corresponding to a radial inflow toward the axis; thus 
(—P), rather than P, diminishes toward the axis. A ra¬ 
dial cross section shows Q increasing outward from the 
axis, until it reaches a maximum at / ^gjris — 1 / 6 ; this 
marks the transition into the envelope region. 

The envelope region follows an additional, geometric 
self-similarity (GSS), such that the two-dimensional flow 
in the ^-77 plane becomes essentially one-dimensional 
when written in terms of the similarity variable x nr 
(defined more precisely in Eq. [34|) . Far from the head. 
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the shock power-law index p = d{ln^s)/d{lnr]) asymp¬ 
totes to yd 2.02; solving the GSS equations (|56H39|) for 
this case (Figure [2j) reproduces the envelope structure. 
Moreover, the P profile, and the external (/ > 3/4; de¬ 
fined by the minimal value of Q in a radial cross section) 
profiles of Q and 1/, are well approximated by the ana¬ 
lytic solution we derive (Eqs. |40] and Figure [T]) for the 
parabolic case yd = 2. 

The total self-similar energy extrapolated from the nu¬ 
merical solution, Etot — 2.1, fixes the dimensionless con¬ 
stant in Eq. (j2T]) as C = 0.8. This is about 

half the value C = 1.5 estimated in GOT for an initial 
jet opening angle of 0.1 radians. The difference is proba¬ 
bly due to the large fraction of the energy deposited far 
beyond the head region. Recovering this energy requires 
filling the attractor and resolving the far downstream, 
which is challenging in a numerical simulation. This in¬ 
complete convergence upon the attractor is likely to be 
the reason for the differences in detail, seen in Eigure[6l 
between our numerical solution and the GOT simulations. 

Our numeric solution appears to be converged, phys¬ 
ical, and unique, because (i)\i qualitatively agrees with 
the simulations of GOT, as provided for the head region 
(Eigure [6]), and with the U sign change reported far 
downstream; (ii)\i agrees qualitatively with the semi- 
analytic constraints of 331 (iii)\i quantitatively agrees 
with the GSS envelope solution, which extends to ^ ^ oo; 
/zr’/convergence tests suggest that it extends to large 
out to our resolution limit (Table [Hand Eigure [5]); and 
(v)iio other solution is found in systematic scans of the 
shock profile at low order {e.g., Eigure |4]), performed in 
two different methods, nor in high order scans. 

Nevertheless, we cannot prove that the solution is an 
attractor, or that it even survives in physical situations 
in which the far downstream transitions into the sub- 
relativistic regime, where our q ^ 1 assumption fails. 
However, the simulations of GOT suggest that a self¬ 
similar solution indeed exists, and is a physically rele¬ 
vant, strong attractor. Under this assumption, our re¬ 
sults substantiate the existence of the solution, uncover 
its structure, in particular far from the head, and show 
that it is unique. 


Moreover, the GOT simulations then imply that the so¬ 
lution is relevant for the typical parameters of external 
GRB shocks, and is probably stable, at least under ax- 
isymmetric perturbations. While computing the obser¬ 
vational implications is straightforward, and a stability 
analysis of the self-similar solution appears feasible, we 
defer these to a forthcoming publication. 

The self-similar regime is strictly valid in the limit of 
an initially extremely narrow and ultra-relativistic jet, 
such that causal contact is reached when the shock is still 
highly relativistic, i.e. the opening angle 0 ^ 1. 

Eor finite initial opening angle Oi and Lorentz factor F^, 
the attractor is only partly filled. The simulations of GOT 
show that explosions with Oi ^ 0.01-0.1 and F^ ^ 20 
indeed produce only a partial attractor, gradually filling 
up throughout the entire quasi-self similar stage. 

For GRB afterglow jets, simulations with Oi > 0.05 


typically show a slow, non-exponential deceleration 


20091: Ivan Eerten et al. 12010 

: Meliani & KeoDensI 

2010: 

van Eerten & MacEadvenI 120111: Wveroda et al.l 

2011; 

De Golle et al.l l2012lj. nnliki 

3 the self-similar evolution 
a self-similar, exponential 
small, Oi < 0.05 initial 

discussed here. However, 
expansion is expected for 

opening anerles (|Wve:oda et al.l 120111: iGranot & PiranI 
1201211 . orovided that the resolution is sufficiently high 

(iGannizzo et al.ll2004 lGranotll200Tl: iMeliani et al.lboOTl): 


such behavior was indeed reported {e.g., Ivan Eert^ 

I2013L and reference therein) for jets with ^ 0.04. 
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APPENDIX 

A. NO REAL-VALUED GSS SOLUTION FOR f3 <2 

Consider the case where oc r]^ far from the axis, with /3 < 2. The GSS scaling is constrained here by the boundary 
conditions (|2S|) . giving rise to modified GSS coordinates Q, P, and U, defined by 

Q = A~'^ri~^Q{x) , P = , U = , and x = C/(^T) > (Al) 

such that Eq. (|2^ becomes 

= m = and U(l) = /3. (A2) 

Far from the axis, the hydrodynamic Eqs. now become 

[1 + 4(/3 - 2)xQ]r + 4(/3 - 2)(2Q + xQ')P = 0 , (A3) 

[1 + (/3 - 2)xQ]Qr + [(/? - 2)Q2 - Q']P = 0 , (A4) 

and 

(U - x/3)r -[P- 2U' + 4(/3 - 2)(U - xU')Q]P = 0 . (A5) 

The first two equations are independent of U, and may be solved analytically; Q is then given by the transcendental 
equation 

[1 - (2 - f3)xQfQ = . (A6) 

For finite values of x, this equation has no real-valued solutions for Q. 

B. A GSS ENVELOPE WITH /3 = 2 GANNOT BE MATGHED TO THE AXIS 

For (3 = 2, the boundary conditions give Us oc r]. In the regime 1/4 < Qao < 1/2, this can be matched to an 

expansion of U along the axis, 

u{i, v) = ciV + + ... , (Bl) 

where are constants that depend only on QaO- In this regime, the expansion leads to 

- -2(2 + QaO) - 

Q(e,^) = r'Q(x) and Qix), (B2) 

where Q is an unknown function that satisfies Q{x cx) = QaO- As x is constant on the shock, Q and P scale 
differently there, contradicting the boundary conditions p6|) . This rules out a monotonic jet with a /3 = 2 GSS 
envelope. 

C. A WIDE JET LEADS TO A NEAR-HEAD DIVERGENGE 

The uniqueness of the numerical solution derived in 0 suggests that the shock profile is fixed by the regularity of 
the flow. To demonstrate how the flow diverges for any deviation from the true shock profile, we consider a wide jet, 
in which an axial expansion series converges rapidly near the head, in the region where the shock is nearly parabolic, 

6 - ih+Ar]'^- 

We expand Qa(0 fo orders a high as n = 5 near the head. For wide jets with A < 0.1, this expansion converges 
rapidly with n near the shock, and shows that Qa becomes non-physically zero and subsequently negative near the 
head. Moreover, at the point where Qa 0, Pa and its derivative must also vanish, and Ui diverges. This vanishing 
of Q near the head of a wide jet is verified numerically, and traced even for A < 0.3. 


D. MONOTONIG HEAD MODEL 


An approximate, monotonic head description may be found using the axial expansion (|43l[45]) . truncated beyond 
order r]^. A monotonic profile far from the head requires Qa(C ^ 1) oc (see ^3.5p . For simplicity, we assume 
that Qa oc even near the head, so the boundary conditions fix Qa = 1/(8^). The axial analysis then implies that 
Pa = ( 1 / 6 )^“^/^. The axial analysis also implies that Ui{^) = 1 , but it is more accurate to determine Ui directly from 
the shock boundary conditions, as shown below. 

Gombining Eqs. (|23H26]) and (|29]), the expansion coefficients of P and Q are fixed up to order r]^, and the coefficients 
of U are fixed up to order 77 ^, by the shock boundary conditions on Q, P and U, their first derivatives, and the second 
derivative of U. 

This can be done for any shock profile, but the resulting coefficient expressions are in general lengthy. For brevity, 
we provide the expansion coefficients for the parabolic profile P accurate close to the nose of the jet. 

Here, the boundary conditions are explicitly given by 


Qs = 


1 

8 ^ 2^2 ^ 8 (^^^2 _|_ _ 8 A 77 ^ 


and 


Ari {A^rf - + 2 ^) 

4 + ihf 


(Dl) 









and 


_ (1 - 2A)A^e + A{A{A{4A - 7) + 7) - 2)^^ - 27l(^ - 1)3^2 
8e(7i(e-a)+a)3 

_^2 ^^3^2 + (^((3 _ 2^)^ _ 5) + 2)^a + (71 - ini) 

8^(^-a)(7i(e-a)+a)3 

.^,71 (-6A3(^ - a)" - 7i2 (^3/3 + isa) (^ - a)" + Tia (19^®/" - isa) (^ - a)) 

i8C3/3(c-a)(7i(c-a) + a)3 

A^h (13^3/3^;, - 7^3/3 - 6C^) 

18f/3(^_^,)(^(^_^,) + ^^)3 ’ 

^ A2 (3A3(C - a)3 + A2 (4^3/3 + 9 Ch) (^ - 4)2 _ ^4 (is^/s _ 94 ) _ 4 )) 

18^3/3(^_4)2(^(^_^4+^43 

A 24 (- 10 ^ 3/34 + 7^3/3 + 3 ^ 1 ) 

+ 184/3(^_4)2(71(^_4)+4)3 ’ 

^2A (A3(e - 4)3 + 15A24(e - 4)2 - A4(4e - 134)(e - a) + 3^^) 

3(A(^-4) + a)3 

r ('■^ =8A3(^ - 4)(A(A^ - {A + 3)4) + 24 ) 

3(A(C-C.)+a)3 

p..^_47l"(A(A + 3)-2)4-4A6^ 

3(A(e-4) + a)3 
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